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Abstract 

In many applications, it is of interest to study trends over time in relationships among 
categorical variables, such as age group, ethnicity, religious afBliation, political party 
and preference for particular policies. At each time point, a sample of individuals pro- 
vide responses to a set of questions, with different individuals sampled at each time. In 
such settings, there tends to be abundant missing data and the variables being mea- 
sured may change over time. At each time point, one obtains a large sparse contingency 
table, with the number of cells often much larger than the number of individuals being 
surveyed. To borrow information across time in modeling large sparse contingency ta- 
bles, we propose a Bayesian autoregressive tensor factorization approach. The proposed 
model relies on a probabilistic Parafac factorization of the joint pmf characterizing the 
categorical data distribution at each time point, with autocorrelation included across 
times. Efficient computational methods are developed relying on MCMC. The methods 
are evaluated through simulation examples and applied to social survey data. 

Key words: Dynamic model; Multivariate categorical data; Nonparametric Bayes; Panel 
data; Parafac; Probabilistic tensor factorization; Stick-breaking. 



1 Introduction 

Time-indexed multivariate categorical data are collected in many areas, with partially- 
overlapping categorical variables measured for different subjects at the different time points. 
As a motivating application, we consider social science surveys that are conducted at reg- 
ular time intervals, containing many categorical questions such as gender, race, age group, 
ethnicity, religious affiliation, political party and preference for particular policies. For such 
surveys and other types of time-indexed multivariate categorical data, it is common for 
the variables measured (questions asked) to vary somewhat over time while a subset of the 
variables will be measured at all times. In addition, the number of variables measured can 
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be moderate to large leading to a contingency table with an enormous number of cells, the 
vast majority of which are empty. Given the fact that social science data often contain com- 
plex interactions, it becomes extremely challenging to build realistic and computationally 
tractable models that allow ultra-sparse data. We define ultra-sparse contingency tables as 
having exponentially or super-exponentially more cells than the sample size. 

Let Xfi = {xtii, ■ ■ ■ jXtipY denote the multivariate response for the ith subject in the 
survey at time t, with the jth categorical question having dj elements, xuj € {1, ... , dj},j = 
1, . . . ,p. We accommodate the case in which the specific variables measured can vary across 
time by introducing missingness indicators, mti = {mtn, . . . , mup)' , with muj = 1 if variable 
j is missing for subject i at time t; we allow design-based missingness in which certain 
variables are not measured for any subjects at a particular time and for individual-specific 
missingness in which certain individuals fail to answer all the questions posed to them. In 
both cases we assume missing at random. 



There is a rich literature on the analysis of contingency tables (jAgrestil (|2002l ) ; 
(j2007l )). Log linear models are perhaps the most commonly used modeling framework. 
Routine implementations rely on maximum likelihood estimation, though there is also a 
rich Bayesian literature. For large, sparse contingency tables, maximum likelihood esti- 
mates do not exist in many cases except for overly-simplistic log-linear models and richer 
classes of models become challenging to implement computationally. There is a rich lit- 
erature on graphical modeling ap proaches to estimating con ditional independence struc- 



Fienberg and Rinaldo 



tures in categorical variables, with 



Dobra and Lenkoski 



( 2011 



proposing a recent Bayesian 



approach. Although their method is computationally efficient, except for very small ta- 
bles, the number of possible graphical models is so enormous that is becomes infeasible 
to visit more than a vanishingly small fraction of the models making accurate rnodel s elec- 



tion or averaging difficult. To fac ilitate scaling to large tables, 



Bhattacharva and DunsonI (^2011 



Dunson and Xind (j2009l ) and 



recently proposed Bayesian probabilistic tensor factoriza- 
tions. These methods express the probability tensor corresponding to the joint probability 
mass function of the categorical variables as a convex combination of independent compo- 
nents. Such methods have not yet been developed for time-indexed contingency tables. 
There is a rich literature on categorical time series and longitudinal data analysis in 
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which the same categorical variable is repeatedly measured for each subject over time. 
For example, Markov models, state space models and random effects models are routinely 
applied in such settings. However, these models are not relevant to the problem of incor- 
porating dependence over time in modeling of large sparse contingency tables. As different 
subjects are measured at different times, we are not faced with the problem of incorporating 
within-subject dependence in repeated observations; instead our goal is to include depen- 
dence in the parameters characterizing the time-dependent joint pmfs for the categorical 
variables. To our knowledge, this problem has not yet been addressed in the literature. Al- 
though one can potentially adapt log-linear or graphical models developed for a contingency 
tables at one time in a somewhat straightforward manner, the hurdles mentioned above for 
the static case are compounded in the dynamic setting. 

To facilitate routine implementations in ultra sparse cases, we propose to adapt the 
Dunson and Xing (DX) (2009) probabilistic Parafac factorization to the dynamic setting. 
The DX model induces a tensor factorization through a Dirichlet process (DP) mixture 
of product multinomial distributions for the categorical observations. There is an increas- 
ingly rich literature proposing nonparametric Bayes dynamic models, which allow time- 
indexed dependent rando m probability measures. Perhaps the most common approach 



relies on a dependent DP (jMacEachernI (|l999l . 



2000)), which incor porates time dep e ndenc e 
in the weights and/o r atom s in a stick-breaking repres entation (jGriffin and Steell (|2006l ): 



Rodriguez and HorstI (|2008l ): IChung and DunsonI (j201ll )). Most applications of dependent 



DPs fix the weights and allow the atoms to vary, as varying weights can lead to computa- 
tional complexities. For dynamic modeling of contingency tables, it is more parsimonious 
to allow varying weights and varying atoms can lead to a substantial computational burden. 
An alternative approach, which allows varying weig hts in a computationally cq i ivenie nt and 



flexible manner, relies of dynamic mixtures of DPs (jDunsonl (1200611: 



cently, a class of probit stick-breaking processes was proposed 



Renetal 



20m)). R e- 



Chung and DunsonI toom 



which has the appealing feature of allowing one to indu ce time dependence in rando m 



probability measures through Gaussian time series models (jRodriguez and DunsonI (|201ll )) 



We propose a new nonparametric state space model for time- indexed ultra sparse contin- 
gency tables. Relying on a DX-type probabilistic Parafac factorization, we place a dynamic 
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model on the weights, which reUes on transformed normal random variables in a similar 
manner to probit stick-breaking. The model is nonparametric in the sense that the in- 
duced prior for each time-indexed joint pmf assigns positive probability in arbitrarily small 
neighborhoods of any "true" data-generating pmf. Hence, our model can allow higher-order 
interactions and complex dependences, while shrinking towards a low-dimensional structure 
and borrowing information across time to address the curse of dimensionality. In addition, 
and crucially for the approach to be useful in the motivating applications, posterior com- 
putation can be implemented via a highly ef ficient Markov ch ain Monte Carlo (MCMC) 



algorithm relying on a slice sampler related to iKalli et al 



(|201ll ). Finally, the factorization 



produces a low-dimensional representation of the joint pmf, which is otherwise characterized 
by a daunting number of parameters in many cases, as the number of cells of the tables can 
be truly massive. 

2 Model specification 

2.1 Modeling of multivariate categorical data 



We review the nonparametric Bayes approach of iDunson and Xind (j2009|) for a static large 
sparse contingency table. Let Xj = {xn, . . . ,Xip)' be multivariate categorical data for the 
zth subject, with Xij E {1, . . . , dj}, j = 1, . . . ,p. Let 

TT = {TTci-Cp, Cj = 1,. . . j = 1, . . . ,p} £ Ud,...dp 



be a probability tensor where 'n'ci-c„ = P{xii = ci, , 



Cp ) is a cell probabil it y and 



^di -dp is the set of all probability tensors of size di x • • • x dp. iDunson and XineJ (|2009l ) 
show that any tt € Hdi - dp can be decomposed as 



TT 



h=l 



where u = (z^i, . . . , v^)' is a probability vector, € Ildi---dp and tp 
dj X 1 probability vec tor for /t = 1 , . . . , fc and 7 = 1, 



(V', 



(i) 

hi ' 



hd-i 



(1) 



IS a 



i2. 



This expression relies on a Parafac 
tensor factorization ()Harshmanl (|l97Cll ) and iKoldal (|200lh ) . It follows that any multivariate 
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categorical data distribution can be expressed as a mixture of product multinomials, 

k p 

P{Xii =Ci,...,Xip = Cp) = TTci-Cp = X] n ^icj- 

h=l j=l 

By introducing a latent class index Si € {1, . . . ,k} for the ith subject, the multivariate 
responses Xj = (xu , . . . , Xjp)' ar e cond itionally independent given Sj. Instead of conditioning 
on a fixed k, Dunson and Xing ( 200^ ) developed a nonparametric Bayes approach that lets 



h=i 



tp^j^^ '-^ Dirichlet(aji, . . . , aja^), independently for j = 1, . . . ,p, 

h = 1, . . . ,oo, 



i^h = Vhllil-Vi), 

l<h 

Vh ~ beta(l, a), independently for = 1, . . . , oo, 



where Cji > for I = 1, . . . ,dj and a > 0. Although (2) allows infinitely many components, 
the number A;„ occupied by the n subjects in the sample will tend to be A;„ << n, so few 
components will be occupied. The model corresponds to a Dirichlet proces s mixture of 



produ ct multinomial distributions relying on a stick-breaking representation (jSethuraman 



(|l994l )). A prior is induced on the joint pmf which has large support in the sense of assigning 



positive probability to Li neighborhoods of any true joint pmf. 

2.2 Modeling of time-indexed multivariate categorical data 

Relying on the DX type probabilistic Parafac factorization, we propose a new nonparametric 
Bayes approach for time-indexed large sparse contingency tables. In a dynamic setting, we 
obtain the time-indexed multivariate response xu = {xui, . . . , xup)' , xuj € {!,..., dj}, for 
the zth subject at time t for i = 1, . . . ,nt, t = 1, . . . ,T and j = 1, . . . ,p. At time t we have 
a probability tensor lit for the multivariate categorical response given by 

T^t = {nci-cp, Cj = l,...,dj, j = l,...,p} G Ud^.-dj, 

where 7Ttc-i---cp = P{xtii = ci, . . . ^xup = Cp) is a cell probability at time t. Relying on the 
probabilistic Parafac factorization, each probability tensor ttj can be expressed as a mixture 
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of product multinomials 

kt 



h=l 



where h e N, Ut = {ua, . . . ,i'tkj is a probability vector, "^th e 11^^...^^, and ip^j^ = 
("^thi^ • • • ' "^thd-y 3l dj xl probability vector for h = I, ... ,kt. Letting su € {I, . . . ,kt} de- 
note a latent class index for the ith subject at time t, the observations x^j are conditionally 
independent given su. 

To borrow information across time, we place a dynamic structure on the probability 
tensor tt^ in ^ assuming time varying weights and static atoms ip^j^ = V'^'^- Time de- 
pendence is induced in the weights through a state space model, which assumes that stick- 
breaking increments on f^/i arise through transforming Gaussian autoregressive processes us- 
ing a monotone differentiable link funct ion q : ^ (0, !)• T h is cha r acterization is motivated 



by th e probit stick-breaking process (IChung and DunsonI (120091 ): [Rodriguez and Dunson 



(j201ll )). and leads to a parsimonious but flexible characterization of time-dependence in 
joint pmfs underlying large, sparse contingence tables. 

Similarly to expression (2), we develop a nonparametric Bayes approach that sets the 
number of components to kt = oo, though the number of occupied components will tend to 
be much less than the sample size and can vary across time. The specific model is 



h=l 

,(i) 



ipf^ ~ Dirichlet(aji, . . . , ajd^), independently for j = 1, . . . ,p, (5) 

/i = 1, . . . , oo, 

iyth = g{Wth)ll{l-g{Wu)}, (6) 

Kh 

Wth = ath + eth, eth^ N{{),al), (7) 
ath = fJ- + (pat-ih + r]th, ~ iV(0, cj^), (8) 

where |(/)| < 1, {eth} and {r]th} are sequences of independently normally distributed random 
variables with mean and variance cr^ and cr^ respectively. The parameter cp controls 
the autocorrelation over time in the weights I'th on the different components. For sake of 
parsimony and simplicity in modeling and computation, we include a single time-stationary 
correlation parameter (p instead of allowing dependence to be time or element specific. In 
the limiting case in which i;^ = 0, the weights uth will be modeled as independent. This 
does not mean that independent priors are placed on the unknown joint pmfs at each 



6 



time, as the incorporation of common atoms automatically induces some degree of a priori 
dependence. However, in applications one typically expects that the joint pmfs will be 
quite similar over time, and by using varying weights one does not rule out arbitrarily 
large changes in the pmfs over time. When ^ is close to one, there will be very high 
time dependence in the weights, leading to effective collapsing on a model that assumes a 
single time stationary joint pmf. For the initial state variables, we assume the stationary 
distributions, aih ^ N{ij,/{1 — 4>),a^/{l — (f?)) independently for ^ = 1, . . . , oo. Also, we 
choose priors /x ~ iV(jLto, o-g), (?!) ~ C/(-l, 1), cr^ ~ IG{me/2, Se/2) and ~ /G(m^/2, 5^/2) 
respectively. 

Expressions (4)- (8) induce a prior on the time-dependent joint pmfs, but it is not im- 
mediately obvious how the chosen hyperpriors in the hierarchical specification impact the 
properties of the prior for {nt}. In particular, it is important to obtain characterizations of 
the moments of the induced prior for the cell probabilities, as well as the prior covariance 
between different elements and across time. Such expressions arc provided in Lemma 1, 
with the proof provided in Appendix A. Lemma 2 shows that the prior is well defined in 
the sense that J^hLi ^th converges to one almost surely. 

Lemma 1. The expectation, variance and covariance of the joint prior on the elements 
of {-Kt} induced through (4)-(8) are 

where /3i = E{g{Wth)}., p2 = ^{^'(H^tfe)}, Ik = E {g{Wth)g{Wt+kh)} , % = E£i %7 and 
!(•) is an indicator function. 

The expectation of cell probabilities can be expressed as the product of expectations 
of Dirichlet priors for atoms. The variance and covariance are expressed as the product of 
two terms, the first one is related to atoms and the second one comes from time varying 
weights. As /i ^ 00, then /?2/(2/3i — 132) 1 and 7fc/(2/5i — 7;^) 1, and the variance 
and covariance will be influenced only by atoms. In such a case, the measure corresponding 
to the stick-breaking process will become a point mass at a random atom almost surely. 
In addition, /5i, (32 and 7^ do not depend on time t, hence all expectation, variance and 
covariance are independent of t though the covariance depends on the time difference k. 
Also, the covariance between cell probabilities with Cj = c'- for all j is always positive and, 
on the other hand, those with Cj ^ for all j have negative covariance. In a special case 
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in which the hyperparameters in the Dirichlet prior are aji = ■ ■ ■ = aj^. = a the variance 
and covariance is zero in the hmit as a — )• oo. The proof is in Appendix A. 
Lemma 2. Yl'h'=i '^th = ^ almost surely. 

Lemma 2 is important in showing that the prior is well defined. The proof is in Appendix 

B. 

Our proposed prior setting is parsimonious but highly flexible in the sense that the 
induced prior assigns positive probability in arbitrarily small neighborhoods of any true 
data-generating pmf. Let H denote the space having elements of the form tt = {tt^ G 
n^i-.-dp, t € {1, . . . ,T}}. We show in Theorem 1 that the proposed prior has large support 

on n. 

Theorem 1. Let Q denote the prior on H through the proposed model and A/'e(7r*') denote 
an Li neighborhood around an arbitrary vr'' € 11. Then for any tt'^ € 11 and e > 0, the prior 
assigns positive probability in the e-neighborhood, Q{j\fe{7^^)} > 0. 

Since the proposed prior is defined o n a space with finitely maii y comp onents, a straight- 
forward extension of theorem 4.3.1 in iGhosh and Ramamoorthil ( 20031 ) ensures that the 
posterior concentrates in arbitrary small neighborhoods of any true data-generating distri- 
bution as the sample size increases. 



3 MCMC algorithm for posterior computation 



For posterior computation in DP mixtures, one common a pproach is marginalizing out t he 
random probability meas ure with the Polya urn sch eme (jBush and MacEachernI (|l996l )). 
Avoiding marginalization, Ishwaran and James ( 200ll ) developed the blocked Gibbs sampler 
relying on truncation app roximation of the stick-breaking rep resentation. Without trun- 



cation, 



Walked (|2007l ) and iPapaspiliopoulos and RobertsI (|2008l ) proposed the slice sampler 



and retrospective MCMC methods respectively. Though the slice sampler is simpler to 
imple ment, conditional constraints on sticks can cause slow mixing of the chain. 



Kalli et al 



(2011 



proposed a more efficient slice sam pler avoiding such a mixing problem. 

(j201ll ). we developed a simple and 



Relying on a slice sampler related to 



Kalli et al 



efficient MCMC algorithm for the proposed model. In the motivating application, we have 
two types of missing data, design-based missingness and individual-specific missingness. We 
assume missing at random for both cases and handle the missing data using missingness 
indicators, mu = {mtn, . . . ^mup)' , with muj = 1 if variable j is missing for subject i at 
time t. In addition, we introduce latent variables ut = {uti, ■ ■ ■ lUtnt)' for the slice sampler. 
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The likelihood of {ut} and {xi} given {rriti}, {vt} and can be expressed as 



T nt 

nn 

t=l i=l 



CO "-J _. / J-. 

j: mtij=0 1=1 



h=l 



This representation is consistent with the original model setting if latent variables {ut} 
are marginalized out. In a special case in which 5 is a probit link function, the data 
augmentation approach in lAlbert and ChibI (j200lh can improve efficiency of the posterior 
sampling by introducing independent normal latent variables {ztif^} with mean Wth and 
variance 1 satisfying 

P{zuh > 0, ztu <0,l<h) = <^{Wth) - HWth)} = yth = P{sti = h). 

l<h 



We propose the following MCMC sampling steps: 

For h = 1, . . . ,k* , with k* = max 
conditional posterior distribution, 



1. For h = 1, . . . ,k*, with k* = maxjstj}, update from the following Dirichlet full 



Dirichlet aji + ^ l{xtij = 1), . . . , aja^ + ^ l{xtij = dj) . 

where Ajh = {{t,i) : muj = 0, su = h}. 
2. Update Zuh from the marginal (w.r.t. uu) conditional posterior distribution. 



ztih • • • ~ < 



N^iWth,!) h<sti, 
N+iWth,l) h = su, 



where N^(Wth,^) and N^(WthA) denote the normal distributions with mean Wth 
and variance 1 truncated on (— cx3,0] and (0, 00) respectively. 

3. Update Wth from the normal marginal (w.r.t. uti) conditional posterior distribution, 
N{Wth,CT'^^J where 



Wth = o-Wth I ^iih + cr^ "^ath I , o-vFi 

A:sti>h 



th sr^nt 



-2 • 

e 



4. Update Uti from the full conditional distribution, Uniform(0, vt 
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5. Update su from the multinomial full conditional distribution, 



Pr{su = h I 



l{h£Bu)U 



-0 V', 



hx 



tij 



^leBti Y[j:mtij=0 '^Ixtij 



where Bu = {h : uth > uu}. To identify the elements in {Ba}, we first update 
ath and Wth for t = 1, . . . , T and h = 1, . . . , fc where k is the smallest number with 
Y!h=i ^th>'^- min{sfi} for all t. 

6. For h = 1, ■ ■ ■ , fe*, update ath usi ng t he forward fjlterixig bac kward sampling algo- 



rithm bv lFriiwirth-Schnatterl (Il994l'l and 



the si mulation smoother by 
(12002). 



Carter and 
de Jong and ShephardI ( 



1995|) and 



<!ohnl (1199411. or Kalman filter and 



Durbin and Koopman 



7. Update /i from the conditional posterior, N{fi^,af^) where fi^ = crf^{a (i + cJq ^uq 



+ a, 



and 



I- ^ ELi Et=2(«tfe - (t>at^ih) + (1 + (A) ELi ^2 
^ A;*{r-l + (l + ,^)/(l -</.)} 



k* {T - 1 + (l + ,^)/(l - ,/.)}■ 



^. Update using the independence MH algorithm in which the proposal distribution 
is constructed relying on the mode and Hessian of the logarithm of the conditional 
posterior densities 7r((/>| • • • ). First, we compute (p which maximizes (or approximately 
maximizes) the conditional posterior density. Then, we generated a candidate from a 
truncated normal distribution TA''(_i cj^), where 



31og7r( 



) 



91og7r( 



9. Update from the conditional distribution, IG{rh^/2, 5'e/2) where rh^ = Tk* + 
and Se = EJ=i EtliiWth - ath? + S^. 

10. Update o"^ from the conditional distribution, IG(m^/2, where = Tk* + 



and 5^ = Y!h=i T,t=2i'=^th - - M-ih)^ + (1 - 0^) ELii^i?* " /^/(l " 4>)V + Sr,- 

In a case in which g is another link function, we update Wth using the independent 
MH algorithm, instead of step 2 and 3 above. We generate a candidate from a normal 
distribution relying on the mode and Hessian of the logarithm of the conditional posterior 
densities of Wth- 



k* 
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4 Simulation study 



In this section, we assess the impact of borrowing of information over time by comparing 
our proposed method to static approaches, such as Dunson and Xing (DX) (2009), apphed 
separately at each time on simulated data. First, we simulate time-indexed contingency 
tables from the model shown in expressions ([l])-® with T = 10, P = 20, dj = 4 for all j, 
^ = 0, = 0.8, o"£ = 0.1 and fi^ = 0.8. At the respective time points we generated 120, 110, 
150, 80, 100, 120, 100, 140, 110 and 150 observations, tiny sample sizes compared with the 
number of cells. For prior distributions, we assumed tpj^^ ~ Dirichlet(l, . . . , 1), /i ~ N[0, 1), 
(j) ~ [/(-1,1), 0-2 ~ /G(2.5, 0.025), ~ /G(2.5, 0.025). We draw 60,000 MCMC samples 
after the initial 20,000 samples are discarded as a burn-in period and every fifth sample 
is saved. We observed that the sample paths were stable and the sample autocorrelations 
dropped smoothly. Therefore, the chains apparently converged and mixed rapidly. 

We first assess performance in estimation of cell probabilities. We picked several cells 
randomly and report true values, posterior means and 95% credible intervals in Figure [J 
(the proposed method) and Figure [2] (DX method). The proposed approach covers all true 
values in 95% intervals and interval widths are much narrower than for the DX approach 
consistently across time. 

We additionally investigate performance in estimating associations among the categor- 
ical variables using the following measure of dependence from Dunson and Xing (2009) 

1 '''' (ncc-i^imPY 

mm{d,,df}-l ' ^ ' 

^ J' J i cj=icj,=i y^tcjVtcj, 

where -ip^p = P{xtij = /) ~ Y^^i=i ^thi^hi ■ The first row of Figure [3] reports plots of all pairs 
of true values (y-axis) and posterior means (x-axis) of ptjj' at time t = 2 and 7. At each 
time point, coordinate points by our approach locate closely to the y = x line, compared to 
widely scattered points by the DX method. In addition. Table [1] shows correlations between 
true values and posterior means of ptjj'- Although correlations by the DX method are high, 
the proposed method consistently produces higher correlations. 





t=l 


t=2 


t=3 


t=4 


t=5 


t=6 


t=7 


t=8 


t=9 


t=10 


Total 


Proposed 
DX 


0.948 
0.837 


0.977 
0.794 


0.990 
0.880 


0.977 
0.761 


0.983 
0.766 


0.986 
0.921 


0.985 
0.846 


0.965 
0.817 


0.969 
0.831 


0.968 
0.793 


0.974 
0.841 



Table 1: Correlations between true values and posterior means of ptjj' using the first simu- 
lation data. 

Log linear models provide a standard choice for the analysis of contingency tables. 
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However, one issue is that flexible log-linear models that accommodate arbitrary interac- 
tions among the variables and allow time dependence cannot be applied directly to large, 
sparse tables. Certainly, maximum likelihood estimates typically do not exist and Bayesian 
methods that a ll ow an unknown dependence structure do not scale beyond small tables. 



Dahinden et al. 



(|201Cl ) proposed an approach for high-dimensional log-linear models with 



interactions, which relies on solvin g several low-di n iensio nal subproblems that are then 
combined. An earlier approach by iDahinden et al.l (120071) instead re li ed on LI penalized 



log-linear models allowing sparsity of tables. Also, iDahinden et al.l (j2007l ) proposed an 
efficient estimation algorithm for model selection for two level categorical variables. 

As a second alternative to our proposed approach, we implemented the method of Dahin- 
den et al. (DH) (2007) in a second simulation example with T = 8, P = 13 and dj = 2 for 
all j. Other settings are the same as in the first simulation case. As DH did not consider 
time-indexed contingency tables, we applied their approach separately at each time point 
using the logilasso R package, with 5-way cross validation used to choose penalty parame- 
ters. The second row of Figure [3] and Table 2 summarize the resulting dependence measures 
Ptjj' at time t = 2 and 7 for each method. For the proposed method, the posterior means 
are close to true values and correlations between estimates and true values are uniformly 
high. The DH method has a tendency to underestimate dependence, particularly when true 
values are low, and has the lowest correlation between the estimates and truth. 





t=l 


t=2 


t=3 


t=4 


t=5 


t=6 


t=7 


t=8 


Total 


Proposed 
DX 
DH 


0.951 
0.872 
0.705 


0.978 
0.803 
0.557 


0.979 
0.838 
0.733 


0.984 
0.599 
0.466 


0.986 
0.807 
0.725 


0.969 
0.884 
0.506 


0.981 
0.932 
0.763 


0.944 
0.827 
0.487 


0.965 
0.696 
0.562 



Table 2: Correlations between true values and posterior means of ptjj' using the second 
simulation data. 

Finally, to gauge robustness we also simulated data from a time-dependent log-linear 
model in which all the coefficients of the main effects and interactions between two variables 
independently follow random walk processes with variance 1 and other higher interactions 
are zero. The third row of Figure [3] and Table [3] report the estimation results. Although 
we find less difference among them in this case, the proposed method still shows the best 
performance. 





t=l 


t=2 


t=3 


t=4 


t=5 


t=6 


t=7 


t=8 


Total 


Proposed 


0.725 


0.827 


0.768 


0.798 


0.818 


0.916 


0.791 


0.807 


0.817 


DX 


0.642 


0.640 


0.726 


0.664 


0.611 


0.864 


0.769 


0.713 


0.724 


DH 


0.371 


0.716 


0.821 


0.491 


0.611 


0.877 


0.764 


0.715 


0.624 



Table 3: Correlations between true values and posterior means of ptjj' using the third 
simulation data. 



12 



5 Analysis of social survey data 



In this section, we apply the proposed method to data from the General Social Sm'vey 
(GSS, http://www3.norc.org /GSS+Website ). Our focus is on studying associations among 



demographic and preference variables over time. We select p = 29 categorical variables from 
1994 to 2010, including gender, ethnicity, preference for particular policies and many more 
listed in the supplemental materials. The GSS was conducted every two years across this 
time period. The numbers of observations are 2,992 (1994), 2,904 (1996), 2,832 (1998), 2,817 
(2000), 2,765 (2002), 2,812 (2004), 4,510 (2006), 2,023 (2008) and 2,044 (2010) respectively. 
There are abundant missing data in which only a subset of the variables were recorded for 
an individual, and compared to the number of cells, the sample size is quite small at each 
time point. 

We first compared our proposed approach to log-linear models. Unfortunately, current 
methodology for fitting log-linear models that allow fiexible dependence structures cannot 
accommodate these data due to the large sparse structure, time variation and abundant 
missing data. Hence, in order to provide a comparison, we initially focused on a bivariate 
subset of the data consisting of religious preference {i = 1,...,5) and attitude towards 
abortion (j = 1, 2) from 1994 to 2010. We consider the following log- linear Poisson models. 

Model 1 : Nuj ~ Poisson (A^t ) , log fuj = A + Af + A/ + A,^^ , 

where Nuj is count of the cell ij at time t, Nt = Nuj, A^ is an effect of the first 

variable (religious preference), A^ is an effect of the second variable (view of abortion) and 
A^^ is an association term. For identifiability, we assume constraints = A^ = A^^ = 
Ag^ = 0. Model 1 assumes no time-dependence in cell probabilities fJ-ij/Y^ii Y^ji ^ii>j>- 

Model 2: Nuj ~ Poisson(iVt fiuj), log ^jLuj = \t + \u + + A^/, 

n — (\ \R \R \A \RA \RA\r 

Pt — [M, \i, ■ ■ ■ , \i > ^tn > • • • 1 \4i ) > 

^ti = f^i + 4>i^t-ii + £th £ti N{0,af), independently for Z = 1, 10, 

where A^, X^j and are effects of the first variable, the second variable and interactions 
at time t respectively. We assume A^ = A^ = A^^ = A^^ = at each time point and 
/3o = for the initial values. Model 2 is a time dependent hierarchical model where all 
parameters in the log-linear model follow AR(1) process independently. 

We firstly estimate all models using the data from 1994 to 2008. Then, relying on the 
estimated parameters, we predict the contingency table in 2010 (Tabled]). For the proposed 
model, we used the same MCMC settings as in the simulation study. For log-linear models. 
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we estimated parameters using an MCMC algorithm where missing values are imputed from 
conditional probabilities given observed data at each iteration. For example, we generate the 
religious preference i given the view of abortion j with probability fiuj/ Yli' fJ'ti'j- For priors, 
we assumed /3 = (A, Af , . . . , Xf, X^, A{\^, . . . , Aff^)' ~ N{0, 1) for Model 1, /i/ ~ iV(0, 1), 
(f>i ~ U{-1,1) and af ~ /G(2.5, 0.025) for aU / for Model 2. Using Gibbs sampling, we 
generated posterior samples of fii and o"| from normal and Inverse-Gamma distributions 
respectively. For (3, (pi, Pt, we used a MH algorithm in which candidates were generated from 
normal distributions relying on the mode and Hessian of the logarithm of the conditional 
posterior densities. We generated 10,000 MCMC samples after the 1,000 burn-in for Model 
1 and 20,000 MCMC samples after the 2,000 burn-in for Model 2 and, for both cases, every 
fifth sample was saved. 

We generated replications at every fifth MCMC iteration and computed average of the 
following predictive criteria. 



5 2 

jobs 
ij 



Absolute deviation (AD): ^ ^ |7V[/^ - 

i=i j=i 

T'rep iKTobs 



1 ^ ^ 

Mean absolute percentage error (MAPE): — 



i=l j = l 

where N^^^ and N^j'"^ are the replication and observation of count of the cell ij respectively. 
To keep the same total number of replications among all methods, predictions are generated 
from cell probabilities fiij/ J2i' Ylj' f^i'j' ^o^' Model 1 and /i20i0jj/ X^j' Ylj' /^20i0i'j' for Model 
2. Table [5] reports the prediction results. Although Model 2 produces better performance 
than Model 1 by incorporating time-dependence, the proposed method clearly outperforms 
log-linear models in terms of both predictive criteria. 





Protestant Catholic 


Jewish 


None Other Total 


Agree 


216 103 


21 


137 60 537 


Disagree 


372 182 


7 


81 47 689 


Total 


588 285 


28 


218 107 1226 


Contingency 


table of the religious 


preference and view of abortion 




Proposed 


Model 1 


Model 2 




AD 194.4 


208.6 


204.5 




MAPE 0.216 


0.232 


0.227 



Table 5: Prediction results. 



Next, we apply the proposed method to all 29 categorical variables. We generated 30,000 
MCMC samples after the initial 10,000 samples are discarded as the burn-in and every fifth 
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sample are saved. We observed the sample paths are stable and the sample autocorrelations 
are small. Table E] shows the estimation result of parameters in the time dependent stick- 
breaking processes. Concerning the measure of time dependence (p, the posterior mean is 
close to 1 and the 95% credible interval locates near 1, which means the weights of the 
stick-breaking processes have strong time dependence over time. 



Parameter 


Mean 


Stdev. 


95% interval 




-0.012 


0.004 


[-0.023, -0.005] 





0.988 


0.004 


[0.978, 0.994] 




0.062 


0.009 


[0.046, 0.082] 


a,, 


0.126 


0.011 


[0.104, 0.149] 



Table 6: Estimation result of parameters in the proposed stick-breaking process. 

Then, we investigate cross interactions among the variables over time. Figure H] show 
the posterior means of ptjj' for all pairs in 2002 and 2010. Additional results for other years 
are included in the supplemental materials. We find the structure of interactions is complex 
at each time point. Also, though each interaction gradually changes over time, all tables 
look similar to one another, implying they have close dependence. This is consistent with 
the result of the strong dependent weights in the stick-breaking processes. Some categorical 
variables such as Race [j = 3], Attitude toward abortion [6], Political party affiliation 
[9] and Think of self as liberal or conservative [14] intricately correlate with many other 
variables. On the other hand, zodiac [11] shows little interactions with all other variables. 
Among all pairs of variables, {Age [1], Marital status [10]}, {Attitude toward abortion [6], 
Attitude toward homosexual [16]} and {Attitude toward homosexual [16], Attitude toward 
Marijuana [19]} show strong interactions in the whole period. Also, we observed several 
pairs of variables showing relatively close interactions over time, such as {Attitude toward 
abortion [6], Think of self as liberal or conservative [14]}, {Race [3], Political party affiliation 
[9]} and {Marital status [10], Having gun [17]}. In addition, the views of government expense 
show moderate interactions, especially to the environment [23], nation's health [24], halting 
the rising crime [25], dealing with drug addiction [26] and education system [27]. 

Next, we study trends of dependence between categorical variables. Figure O reports the 
posterior means and 95% credible intervals of ptjj' for pairs with close interactions. We ob- 
served various patterns of time paths. For {Age, Marital status}, the interaction increased 
around 2000 then declined sharply to a lower level. {Race, Political party affiliation} and 
{Race, Having gun} have peaks in 2006 and the interactions have steeply decreased after 
that. In addition, we can see similar trends in {Attitude toward abortion. Think of self 
as lib or con}, {Attitude toward abortion, Attitude toward homosexual}, {Attitude to- 
ward homosexual, Attitude toward Marijuana}, {Religion, Attitude toward abortion} and 
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{Religion, Attitude toward Marijuana}. The interactions have roughly increased over time, 
especially in the 2000s. On the other hand, the dependence in {Race, Death penalty for 
murder} decreased at first and kept stable in the middle of the period then declined again. 
{Having gun, Family income} gradually increased over the period but the difference is small. 
For {Marital status. Having gun}, the interaction dropped in the middle of the period but 
recovered recently at the same level as the beginning. 

6 Discussion 

We have demonstrated that the proposed approach is useful in analyzing time-indexed large 
sparse contingency tables. One interesting extension is to accommodate joint modeling of 
mixed scale variables consisting of not only categorical data but also continuoTis and count 
variables. In such a case, one can potentially model the observed data vector for the ith 
subject at time t, yu = {uui, ■ ■ ■ lyupY , as conditionally independent given latent class 
variables xu = {xm, . . . , Xup)' , with xu modeled exactly as proposed in this article. For 
example, consider the simple case in which p = 2 with ytn € 3? continuous and yti2 € 
{!,..., ^2} categorical. Then, one can let ym ~ ^{f^xui^^xtn) ^^'^ Vti^ — ^ti2, with the 
proposed probabilistic tensor factorization approach flexibly accommodating dependence in 
ytii and yti2 through dependence in Xfn and Xti2- The induced marginal distribution for 
the continuous variable yui will be a mixture of normals, with the probability weight on 
each component potentially varying with the categorical variable yti2- This same strategy 
can be generalized to more complex settings involving many categorical, count, continuous 
and even functional observations. 

Another interesting direction in terms of generalizations is to accommodate dependence 
in the observations; for example, one may collected multivariate categorical longitudinal 
data in which the same variables are measured repeatedly on the sample study subjects 
or the data may have a nested structure. Log linear and logistic regression- type models 
can be easily generalized to such settings, but clearly encounter computational challenges 
in large sparse settings. Potentially the simplex factor model of Bhattacharya and Dun- 
son (2012) can be generalized to accommodate such dependence structures through the 
latent factors, with some challenges arising in terms of developing computationally efficient 
implementations and models that are both flexible and interpretable. 
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A Proof of Lemma 1 



The expectation of cell probability is 



oo p 



h=i j=i 



E 

h=l 



P CXI p p 

n^{<}E^w=n^{<}=n 



h=i 



The marginal distribution of Wth can be expressed as N{^/{1 — (/)),cr^/(l — <p'^) + cr^), 
independent of t and h. Hence, we set /3i = E{g{Wth)} and (32 = E{g'^{Wth)}- The second 
moment of cell probability is 



^{<...cj = E 



oo p I I oo p 



h=i j=i 



oo oo 



1=1 j=l 



h=l 1=1 

p 



oo p oo oo 



h=l 1=1 



7 = 1 ""J^"^ + i = l '^i / h=l j = l ~o 



'■1 1 

a. 



where 



E^{-'^} = E^ 



h=i 



h=l 
oo 



9\Wth)ll{l-9{Wu)y 



Kh 



J]/32{l-2/3i+/32}"-\ 



h=i 



2/3i - /32 
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Hence, 



P a? 



(10) 



Similarly, 

-E^{7rtci...cp7ri+fcc',-c^} = E 



00 p r 00 p 

/i=i i=i I 1 1=1 i=l 



i=i j=i 

^ flic,- {ojc' + 1 (cj = ) } ^ Ojc, a^c; 



00 p 



h.=l 



n 



aj{aj + 1) 



n 



i=l J / h=l 



" 2 ' 

• 1 



where 



Ell^thJ^t+kh} = E 



9{Wth)\{{l - 9{Wu)} 



Kh 



g{Wt+kh)\{{l - 9{Wt+ki)} 



Kh 



= E {g{Wth)9{Wt+kh)} \{E[{1- g{Wu)}{l - g{Wt+ki)}] , 

Kh 

= E{g{Wth)9{Wt+kh)} n [1 - + E{g{Wti)g{Wt^ui)}] ■ 

Kh 

From ([7]) and (jH]), E {g{Wth)9{Wt+kh)} can be expressed as 

E {9{Wth)g{Wt+kh)} = E {g{ath + eth)g{at+kh + £t+kh)} , 

= £' {ath + eth) g ( ^ _ ^ /U + 4>^ath + ^ <p'wt+k-ih + ^i+fch 1 | • 

Since ajh, £th, wt+k-ih (i = 0, . . . , A; — 1) and et+kh are independent of one another and 
their distributions do not depend on t or h, hence 7^ = E {g{Wth)giWt+kh)} is dependent 
on time difference k but independent of time t. 
In addition, 

00 00 

^ E{uthiyt+kh} = Y,^kll{l -2(3i+ 7fc} , 

h=l h=l Kh 

= 

Wl-lk 
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Hence, 



p 



Since /32/(2/3i — /32) > 0, 7,fc/(2/3i — 7^) > and pUj) . cell probabilities with Cj = c^- for all j 
have positive covariance and, on the other hand, those with Cj ^ c'j for all j have negative 
covariance. 

In a case where aji = ■ ■ ■ = ajcj = a, the variance and covariance are expressed as 
Hence, l^{vrtei...cp} and Co7;{7rtei...cp, as a ^ 00. 

B Proof of Lemma 2 



b pro ve X^^i ft/i = 1 a.s., it is enough to show Yl'h'=i E{^^&{^~9(^th)} = —00 ([ishwaran and James 



(|200ll )). 51 is a non-negative monotone increasing link function: 3? — ?> (0,1), therefore 



< /3i = E{g(Wth)} < 1- Then, using Jensen's inequality, 

^[log{l - g{Wth)}] < log[l - E{g{Wth)}] = log(l - < 0. 
Therefore, Yl'h'=i ^{log(l ~ di^th)} = —00 at each time point. 

C Proof of theorem 

The proposed prior probability assig ned to Me{n°) can be expressed as 
Q{A4(7rO)} = j l{\\7v-7v^\\<e)dQ{iyt,'ip'il\te{l,...,T},h = l,...,^,j = l,...,p). 
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where Vt is a probability vector induced by the proposed stick breaking process and we use 
the Li distance 



TV — TV 



,0 



tC\---Cp I ' 



t=l Cl = l Cp = l 



where is a probabihty mass function for time t G {1, . . . ,T}. 



For any tv € H, each component in tv can be expressed as 



th 1 



h=l 



where /cj G N, i^f = (t'li, • • • , i^f^)' is a probabihty vector, ^'^/j G n^i-dp and ■j/j^^^'^ = 
('^thi^ • • • ^'^^hd )' a X 1 probabihty vector. We define = Q and /c^ = X]i=i ^« ^^^^ ^ = 
1, . . . ,T. Then, we construct tt = {vr^jt G {!,••• G H induced by the proposed prior 

such that the component with the index h in ttJ* is approximated by the component with the 
index kt^ + /i in ttj. Let Ut = (i^ti, i^t2T ■ ■)' be a probabihty vector where Dtm = i^? .+ 
for k^_^ < m < k^ and = otherwise, i.e., = u^f^ for 1 < h < kt- For any 

e, we define a set D{Tv^,e) C 11 such that for any tv G L'(7r'^,e), each ttj can be expressed 
as dH) satisfying G Me'{i>), where = G {l,...,r}}, i> = {i>t,t G {1,...,T}} 

and e' = e/2nj=i c^j, and ^[^^ G K" (^JP) for /i = 1, . . . , fct and t = 1, . . . , T where 
e" = e/2ZtPthpUjdr 

We consider the intervals {ath,bth) in the real line for Wth in the proposed prior for 
h = 1, . . . ,kf and t = 1, . . . , T where 

[5-Hmax(% -e,0)}, (/i = 1), ig'^ihh + e}, {h = l), 



where e = ^' /^Y^f-Ptk^ . In this case, it is straightforward to check lut^ — i'th\ < e for 
h = 1, . . . ,kf and the proposed prior assigns positive probability to these intervals. Then, 
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the distance between v and v is 



T oo 
i=l h=\ 

= ^pt^Wth - i>th\ + J2p^ ^ ^^'^ 

i=l h=l t=l h>k+ 

T 

<2eY,Ptkt = e'. 



(11) 



t=i 



For the second component in (fTTj) . i^t/i < kfe because vth > t'th, — e for /i = 1, . . . , /c^ 



and X^/jLi ^th > 1 — e. In addition, it is straightforward to show that the proposed prior 
assigns positive probabihty to M^" (^fplh^- Therefore, since D{7v'^,e) contains such case, 
Q{D{7vO,e)}>0. 

For any tt € D{7t^, e). 



T di 



t=l Cl = l Cp = l 



T di dp 

t=l Cl=l Cp = l 

T di dp 

t=l Cl=l Cp = l 



oo p kt P 

h=l i=l i=l i=l 



kt 



E|-.e...n''il 



■°nw/.u+ E ""11*. 

l<k+^^,kf<l i=i 



-he, 



/l=l 



fiE«E-E E 

t=l Cl=l Cp = l 



T di 



h=l 

kt 



(i) 



1=1 j=i 



4=1 ci=l Cp=l 
dl dp T oo 



\ 



h=l 



T kt P di 



ci=l Cp=l t=l h=l t=l 1=1 j=l ci=l Cp=l 

p T p 

< Yldje' + Y,Ptktplldje", 

j=i t=i j=i 

e e 

= 2 + 2=^- 

Therefore tt G A/;(7r°) and D{TT^,e) C A4(7r°). Hence, Q{M,{tt^)} > 0. 
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Figures 
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2 4 6 8 10 f 2 4 6 8 10 f 

The first row: P{xtii = 0,Xti6 = i,Xtiio = 2, xtas = 3) and P{xtn = 2,xti9 = 

0,XU13 = 3,XU19 = 1). 

The second row: P{xtii = 2,xu7 = l,a;«2o = 3) and P{xuz = 3,a;«i2 = ^,xtns = 0). 
The third row: P(xtiii = 1, xun = 1) and P{xu5 = 2, xuig = 1). 
The forth row: P{xti» = 0) and P{xti2Q ~ 3). 

Figure 1: Estimation results of cell probabilities by the proposed method. 
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The first row: P{xu4, = 0,Xti6 = i,Xtiio = 2, xtas = 3) and P{xtn = 2,xu9 = 

0,XU13 = 3,XU19 = 1). 

The second row: P{xtii = 2,xm = l,xu2o = 3) and P{xuz = 3,a;«i2 = ^,xtns, = 0). 
The third row: P{xuii = 1, xun = 1) and P{xu5 = 2, xuig = 1). 
The forth row: P{xtis = 0) and P{xti20 = 3). 



Figure 2: Estimation results of cell probabilities by DX method. 
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Case1(t = 2) 




0.05 0.1 0,15 0.2 0.25 0.3 
True values 



Case1(t = 7) 




0.05 0.1 0.15 0.2 0.25 0.3 
True values 



Case 2 (t = 2) 




0.05 0,1 015 0.2 0,25 0.3 0.35 0.4 
True values 



Case 2 (t = 7) 




0.05 0,1 0.15 0.2 0.25 0.3 0.35 0.4 0,45 0.5 
True values 



Case 3 (t = 2) 




0,05 0.1 0.15 0.2 0.25 0.3 0,35 04 0.45 
True values 



Case 3 (t = 7) 




0,1 02 0,3 0.4 0,5 
True values 



y cixis represents estimated values and x ajds true values. Cross-shaped dots represent the proposed 
method, circles DX method and triangles DH method. The first, second and third rows show the results 
at time t = 2 and 7 using the first (case f), second (case 2), third (case 3) simulation data sets. 

Figure 3: Plots of true and estimated values of ptjj' using the simulation data. 
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The first row: (Age group, Current marital status), (Attitude toward abourtion. Think of self as liberal or 
conservative) and (Attitude toward abourtion. Attitude toward homosexual sex relations) . 

The second row: (Attitude toward homosexual sex relations. Should Marijuana be made legal), (Race, Political 
party afRliation) and (Race, Favor or oppose death penalty for murder). 

The third row: (Race, Have gun in home), (Religious preference. Attitude toward abourtion) and (Have gun in 

home, Total family income). 

The fourth row: (Current marital status. Have gun in home) and (Religious preference. Should Marijuana be 
made legal). 



Figure 5: Estimation results of ptjj' for several pairs. 
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E Supplemental materials 



No. Categorical variable (Name in GSS) # of categories 

"1 Age group* (AGE) 8~ 

2 Sex (SEX) 2 

3 Race (RACE) 3 

4 Religious preference** (RELIG) 5 

5 Region (REGION) 9 

6 Attitude toward abortion (ABANY) 2 

7 Should Govetnment help pay for medical care? (HELP SICK) 5 

8 Highest degree (DEGREE) 5 

9 Political party affiliation (PARTYID) 8 

10 Current marital status (MARITAL) 5 

11 Astrological sign (ZODIAC) 12 

12 Confidence in banks and financial institutions (CONFINAN) 3 

13 Confidence in U.S. Supreme Court (CONJUDGE) 3 

14 Think of self as liberal or conservative (POLVIEWS) 7 

15 Belief in fife after death (POSTLIFE) 2 

16 Attitude toward homosexual sex relations (HOMOSEX) 5 

17 Have gun in home (OWNGUN) 2 

18 Subjective class identification (CLASS) 4 

19 Should Marijuana be made legal (GRASS) 2 

20 Total family income (INCOME) 12 

21 Favor or oppose; death penalty for murder (CAPPUN) 2 

22 Attitude toward spending money on space exploration program (NATSPAC) 3 

23 Attitude toward spending money on improving and protecting environment (NATENVIR) 3 

24 Attitude toward spending money on improving and protecting the nations's health (NATHE3\.L) 

25 Attitude toward spending money on halting the rising crime rate (NATCRIME) 3 

26 Attitude toward spending money on dealing with drug addiction (NATDRUG) 3 

27 Attitude toward spending money on improving the nation's education system (NATEDUC) 3 

28 Attitude toward spending money on the military, armaments and defense (NATARMS) 3 

29 Attitude toward spending money on foreigh aid (NATAID) 3 



*The category of Age group is different from the original one: 1. 18 or 19 years old, 2. 20s, 3. 30s, 4. 40s, 
5. 50s, 6. 60s, 7. 70s, 8. more than 80 years old. 

**The category of Religious preference is different from the original one: 1. Protestant, 2. Catholic, 3. 
Jewish, 4. None, 5. Others. 

Table 7: List of categorical variables. 
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Figure 7: Posterior means of ptjj' in 1998 (above) and 2000 (below). 
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